rm(list=ls());gc()
##          used (Mb) gc trigger  (Mb) max used  (Mb)
## Ncells 463965 24.8    1000164  53.5   650900  34.8
## Vcells 852696  6.6   17698584 135.1 20029491 152.9
unlink('Causal_Impact2_cons_trauma_cache', recursive = TRUE)

load(paste0(getwd(),"/","Procesos hasta 4_2.RData"))
#xaringan::inf_mr()

if(isTRUE(getOption('knitr.in.progress'))==T){
    clus_iter=30000
} else {
  input <- readline('¿Are you gonna run the dataset with the whole iterations? (Si/No): ')
  if(input=="Si"){
    clus_iter=10000
  } else {
    clus_iter=1000
  }
}
#arriba puse algunas opciones para que por defecto escondiera el código
#también cargue algunos estilo .css para que el texto me apareciera justificado, entre otras cosas.
local({r <- getOption("repos")
       r["CRAN"] <- "http://cran.r-project.org" 
       options(repos=r)
})

`%>%` <- magrittr::`%>%`
copy_names <- function(x,row.names=FALSE,col.names=TRUE,dec=",",...) {
  if(class(ungroup(x))[1]=="tbl_df"){
        if(options()$OutDec=="."){
            options(OutDec = dec)
            write.table(format(data.frame(x)),"clipboard",sep="\t",row.names=FALSE,col.names=col.names,...)
            options(OutDec = ".")
          return(x)
        } else {
            options(OutDec = ",")
            write.table(format(data.frame(x)),"clipboard",sep="\t",row.names=FALSE,col.names=col.names,...)
            options(OutDec = ",")
          return(x)    
        }
  } else {
        if(options()$OutDec=="."){
            options(OutDec = dec)
            write.table(format(x),"clipboard",sep="\t",row.names=FALSE,col.names=col.names,...)
            options(OutDec = ".")
          return(x)
        } else {
            options(OutDec = ",")
            write.table(format(x),"clipboard",sep="\t",row.names=FALSE,col.names=col.names,...)
            options(OutDec = ",")
          return(x)       
  }
 }
}  

if(!require(pacman)){install.packages("pacman")}

pacman::p_unlock(lib.loc = .libPaths()) #para no tener problemas reinstalando paquetes
knitr::opts_chunk$set(
    echo = TRUE,
    message = FALSE,
    warning = FALSE
)
#dejo los paquetes estadísticos que voy a utilizar

if(!require(plotly)){install.packages("plotly")}
if(!require(lubridate)){install.packages("lubridate")}
if(!require(htmlwidgets)){install.packages("htmlwidgets")}
if(!require(tidyverse)){install.packages("tidyverse")}
if(!require(gganimate)){install.packages("gganimate")}
if(!require(readr)){install.packages("readr")}
if(!require(stringr)){install.packages("stringr")}
if(!require(data.table)){install.packages("data.table")}
if(!require(DT)){install.packages("DT")}
if(!require(ggplot2)){install.packages("ggplot2")}
if(!require(lattice)){install.packages("lattice")}
if(!require(forecast)){install.packages("forecast")}
if(!require(zoo)){install.packages("zoo")}
if(!require(panelView)){install.packages("panelView")}
if(!require(janitor)){install.packages("janitor")}
if(!require(rjson)){install.packages("rjson")}
if(!require(estimatr)){install.packages("estimatr")} 
if(!require(CausalImpact)){install.packages("CausalImpact")}
if(!require(textreg)){install.packages("textreg")}
if(!require(sjPlot)){install.packages("sjPlot")}
if(!require(foreign)){install.packages("foreign")}
if(!require(tsModel)){install.packages("tsModel")}
if(!require(lmtest)){install.packages("lmtest")}
if(!require(Epi)){install.packages("Epi")}
if(!require(splines)){install.packages("splines")}
if(!require(vcd)){install.packages("vcd")}
if(!require(astsa)){install.packages("astsa")}
if(!require(forecast)){install.packages("forecast")}
if(!require(MASS)){install.packages("MASS")}
if(!require(ggsci)){install.packages("ggsci")}
if(!require(Hmisc)){install.packages("Hmisc")}
if(!require(compareGroups)){install.packages("compareGroups")}
if(!require(dplyr)){install.packages("dplyr")}
if(!require(ggforce)){install.packages("ggforce")}
if(!require(imputeTS)){install.packages("imputeTS")}
if(!require(doParallel)){install.packages("doParallel")}
if(!require(SCtools)){install.packages("SCtools")}
if(!require(MSCMT)){install.packages("MSCMT")}
# Calculate the number of cores
no_cores <- detectCores() - 1
cl<-makeCluster(no_cores)
registerDoParallel(cl)

Sys.setlocale(category = "LC_ALL", locale = "english")
## [1] "LC_COLLATE=English_United States.1252;LC_CTYPE=English_United States.1252;LC_MONETARY=English_United States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252"
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#
######1. Ejemplo 1 #####
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#

#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_

#PARÁMETROS INGRESADOS A LOS MODELOS POR THOMAS
#cons_trauma ~ offset(log(offset)) + year + as.Date(date) + month + day + weekday + yearday + prevtrc + difftrc + hosp_trauma,
#hosp_trauma ~ offset(log(offset)) + year + as.Date(date) + month + day + weekday + yearday + prevtrh + difftrh + cons_trauma,
#data15a64[data15a64$tx == 0 & data15a64$txtime == 0, ]
#glm(hosp_trauma ~ hosp_circ + year + month + day + weekday + yearday + prevtrh, family = quasipoisson(), data15a64[data15a64$tx == 0 & data15a64$txtime == 1,])
#glm(hosp_trauma ~ hosp_circ + year + month + day + weekday + yearday + prevtrh, family = poisson(), data15a64[data15a64$tx == 0 & data15a64$txtime == 1,])
#lm(hosp_trauma ~ hosp_circ + year + month + day + weekday + yearday + prevtrh, data15a64[data15a64$tx == 0 & data15a64$txtime == 1,])

#Genero una columna conteo por fila
data15a64_rn<- data15a64_wk%>% dplyr::mutate(rn=row_number())
data65a_rn<-  data15a64_wk%>% dplyr::mutate(rn=row_number())

#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#

rn_summ_data15a64<-
data15a64_rn%>%
    dplyr::filter(did==1)%>%
    dplyr::summarise(min=min(rn),max=max(rn),diff=max-min)
post.period <- c(as.numeric(rn_summ_data15a64["min"]), as.numeric(rn_summ_data15a64["max"]))


1 Bayesian Structural Time-Series


In order to reproduce adequately what was solicited, we looked over the models Thomas did. The independent variables we selected were: - Circulatory Consultations (cons_circ) - Respiratory Consultations (cons_resp) - Year, month, day, weekday, yearday - Difference Between Total Hospitalizations and Trauma Hospitalizations (difftrh)


The model of Bayesian Structural Time Series (BSTS) with Causal Impact Analysis contains a local level component (observation equation linking observed data to a state vector) and a seasonal component (state equation that describes how the state vector evolves over time). Also it can contain other predictor series. The method uses a Markov chain Monte Carlo (MCMC) sampling of the structural time series given the observed data and the priors.


1.1 Several models contrasted

We decided to contrast several models with different specifications and control variables. We started with a model that did not have a counterfactual.


1.1.1 No counterfactual, no control variables

library(CausalImpact)

data15a64_rn2<-
data15a64_rn%>%
  dplyr::mutate(log_offset=log(offset))%>%
  dplyr::mutate(log_cons_trauma=log(cons_trauma))%>%
  dplyr::mutate(log_cons_circ=log(cons_circ))%>%
  dplyr::mutate(log_cons_resp=log(cons_resp))%>%
  dplyr::mutate(log_difftrh=log(difftrh))%>%
  data.frame()
  
#rn_summ_data15a64

#CREATE THE DB --> SUPER IMPORTANT TO ADD THE OUTPUT OF INTEREST AS THE FIRST COLUMN here, plus to select all the other time series we want:
#in this case we take all the columns from 2 to 10
df <- zoo(data15a64_rn2[c("cons_trauma","cons_circ","cons_resp","difftrh")], as.Date(data15a64_rn2$date)) #prevtrh had a missing value due to the absence of previous value
#set the period before intervention
pre.period <- as.Date(c("2015-01-01", "2019-10-17"))
#set the period after intervention
post.period <- as.Date(c("2019-10-18", "2019-12-31"))
#compute CausalImpact, look at the documentation to change niter and nseasons
set.seed(2125)
impact <- CausalImpact(df, pre.period,post.period, model.args=list(niter=clus_iter, nseasons=7))
plot(impact, "original") 

plot(impact)

##summary(impact) 
##summary(impact,"report")


1.1.2 With counterfactuals

We decided to exclude the post-treatment period data to specify our model, to mimic the fact that this data should be unobserved after the intervention. We introduced yearly and monthly seasonal components into the time-series structure. We started with a Random Walk. This local level model assumes the trend is a random walk (do not assume an observable pattern or trend).


###RANDOM WALK
#######https://magoosh.com/statistics/what-is-random-walk-theory/
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_

##data%>% dplyr::filter(yearday>0)%>% group_by(year)%>% summarise(min(as.Date(date)),max(as.Date(date)),n())

#I will explore their cumulative 1 step forward prediction error for the time series between the start to the time of intervention
# Post Intervention Period is filled with NA
# Remove outcomes from the post-period. The BSTS model should be ignorant of the values we intend to predict

#_#_#_#_#_#_#_#_#_#_#_#_#_
#MADE WITH VECTORS: https://rpubs.com/irJERAD/Causal-Impact
#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#c("log_cons_resp","log_cons_circ","log_cons_trauma","log_diffreh")
data15a64_rn_causal <- data15a64_rn %>%
    dplyr::mutate(cons_trauma = replace(cons_trauma, did >= 1, NA)) %>% 
    dplyr::mutate(log_cons_trauma=log(cons_trauma))%>%
    dplyr::mutate(log_cons_circ=log(cons_circ))%>%
    dplyr::mutate(log_cons_resp=log(cons_resp))%>%
    dplyr::mutate(log_difftrh=log(difftrh))

post_period_response <- as.numeric(unlist(data15a64_rn[which(data15a64_rn$did==1),"cons_trauma"]))


y<- data15a64_rn_causal$cons_trauma
x1 <-data15a64_rn_causal$cons_circ
x2 <-data15a64_rn_causal$cons_resp
x3 <-data15a64_rn_causal$difftrh
x4 <-data15a64_rn_causal$log_offset

#y<- data15a64_rn_causal$log_cons_trauma
#x1 <-data15a64_rn_causal$log_cons_circ
#x2 <-data15a64_rn_causal$log_cons_resp
#x3 <-data15a64_rn_causal$log_difftrh
#x4 <-data15a64_rn_causal$log_offset

# Model 1
ss <- list()
# Local trend, weekly-seasonal #https://qastack.mx/stats/209426/predictions-from-bsts-model-in-r-are-failing-completely - PUSE UN GENERALIZED LOCAL TREND
ss <- AddLocalLevel(ss, data15a64_rn_causal$cons_trauma) #
# Add weekly seasonal
ss <- AddSeasonal(ss, data15a64_rn_causal$cons_trauma, nseasons=5, season.duration = 52) #weeks OJO, ESTOS NO SON WEEKS VERDADEROS. PORQUE TENGO MAS DE EUN AÑO
ss <- AddSeasonal(ss, data15a64_rn_causal$cons_trauma, nseasons = 12, season.duration = 4) #months
# For example, to add a day-of-week component to data with daily granularity, use model.args = list(nseasons = 7, season.duration = 1). To add a day-of-week component to data with hourly granularity, set model.args = list(nseasons = 7, season.duration = 24).
model1 <- bsts(data15a64_rn_causal$cons_trauma, 
               state.specification = ss, #A list with elements created by AddLocalLinearTrend, AddSeasonal, and similar functions for adding components of state. See the help page for state.specification.
               #family ="poisson", #A Bayesian Analysis of Time-Series Event Count Data
               niter = clus_iter, 
               #burn = 200, #http://finzi.psych.upenn.edu/library/bsts/html/SuggestBurn.html Suggest the size of an MCMC burn in sample as a proportion of the total run.
               seed= 2125)
## =-=-=-=-= Iteration 0 Tue Feb 09 00:05:31 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 3000 Tue Feb 09 00:05:44 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 6000 Tue Feb 09 00:05:56 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 9000 Tue Feb 09 00:06:09 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 12000 Tue Feb 09 00:06:21 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 15000 Tue Feb 09 00:06:34 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 18000 Tue Feb 09 00:06:46 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 21000 Tue Feb 09 00:06:59 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 24000 Tue Feb 09 00:07:11 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 27000 Tue Feb 09 00:07:24 2021
##  =-=-=-=-=
plot(model1, main = "Model 1")

plot(model1, "components")

impact2 <- CausalImpact(bsts.model = model1,
                       post.period.response = post_period_response)
plot(impact2, "original") 

burn1 <- SuggestBurn(0.1, model1)

plot(impact2)

#summary(impact2) 
##summary(impact2,"report")
#msts(data15a64_rn$hosp_trauma, seasonal.periods=c(7,365.25))
#The response variable (i.e., the first column in data) may contain missing values (NA), but covariates (all other columns in data) may not. If one of your covariates contains missing values, consider imputing (i.e., estimating) the missing values; if this is not feasible, leave the regressor out.

# Model 2
ss2 <- list()
# Local trend, weekly-seasonal #https://qastack.mx/stats/209426/predictions-from-bsts-model-in-r-are-failing-completely - PUSE UN GENERALIZED LOCAL TREND
ss2 <- AddLocalLevel(ss2, data15a64_rn_causal$cons_trauma) #
# Add weekly seasonal
ss2 <- AddSeasonal(ss2, data15a64_rn_causal$cons_trauma, nseasons=5, season.duration = 52) #weeks OJO, ESTOS NO SON WEEKS VERDADEROS. PORQUE TENGO MAS DE EUN AÑO
ss2 <- AddSeasonal(ss2, data15a64_rn_causal$cons_trauma, nseasons = 12, season.duration =4) #years
#ss2 <- AddAutoAr(ss2, y = data15a64_rn_causal$log_hosp_trauma, lags = 1) #NO PUEDO AREGAR AR1 CON POISSON
# For example, to add a day-of-week component to data with daily granularity, use model.args = list(nseasons = 7, season.duration = 1). To add a day-of-week component to data with hourly granularity, set model.args = list(nseasons = 7, season.duration = 24).
model2 <- bsts(data15a64_rn_causal$cons_trauma ~ x1, 
               state.specification = ss2, #A list with elements created by AddLocalLinearTrend, AddSeasonal, and similar functions for adding components of state. See the help page for state.specification.
               #family ="poisson", #A Bayesian Analysis of Time-Series Event Count Data
               niter = clus_iter, 
              # burn = 500, #http://finzi.psych.upenn.edu/library/bsts/html/SuggestBurn.html Suggest the size of an MCMC burn in sample as a proportion of the total run.
               seed= 2125)
## =-=-=-=-= Iteration 0 Tue Feb 09 00:08:17 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 3000 Tue Feb 09 00:08:30 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 6000 Tue Feb 09 00:08:42 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 9000 Tue Feb 09 00:08:54 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 12000 Tue Feb 09 00:09:06 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 15000 Tue Feb 09 00:09:19 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 18000 Tue Feb 09 00:09:31 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 21000 Tue Feb 09 00:09:43 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 24000 Tue Feb 09 00:09:55 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 27000 Tue Feb 09 00:10:08 2021
##  =-=-=-=-=
plot(model2, main = "Model 2")

plot(model2, "components")

impact3 <- CausalImpact(bsts.model = model2,
                       post.period.response = post_period_response)
plot(impact3, "original") 

burn2 <- SuggestBurn(0.1, model2)

plot(impact3)

#summary(impact3) 
##summary(impact3,"report")
#dynamic.regression Whether to include time-varying regression coefficients. In combination with a time-varying local trend or even a time-varying local level, this often leads to overspecification, in which case a static regression is safer. Defaults to FALSE.


1.1.3 Local Linear Trend

A model with too many components can sometimes offer too much flexibility, providing unrealistically widening forecasts. This is why the default model does not include a local linear trend component (see this link). But in the next two models we assumed a local structure of the latent state variable (AddLocalLinearTrend), instead of a random process.


# Model 1
ssb <- list()
# Local trend, weekly-seasonal #https://qastack.mx/stats/209426/predictions-from-bsts-model-in-r-are-failing-completely - PUSE UN GENERALIZED LOCAL TREND
ssb <- AddLocalLinearTrend(ssb, data15a64_rn_causal$cons_trauma) #AddSemilocalLinearTrend #AddLocalLevel
# Add weekly seasonal
ssb <- AddSeasonal(ssb, data15a64_rn_causal$cons_trauma, nseasons=5, season.duration = 52) #weeks OJO, ESTOS NO SON WEEKS VERDADEROS. PORQUE TENGO MAS DE EUN AÑO
ssb <- AddSeasonal(ssb, data15a64_rn_causal$cons_trauma, nseasons = 12, season.duration =4) #years

# For example, to add a day-of-week component to data with daily granularity, use model.args = list(nseasons = 7, season.duration = 1). To add a day-of-week component to data with hourly granularity, set model.args = list(nseasons = 7, season.duration = 24).
model1b <- bsts(data15a64_rn_causal$cons_trauma, 
               state.specification = ssb, #A list with elements created by AddLocalLinearTrend, AddSeasonal, and similar functions for adding components of state. See the help page for state.specification.
               #family ="poisson", #A Bayesian Analysis of Time-Series Event Count Data. NO SE PUEDE OCUPAR
               niter = clus_iter, 
               #burn = 200, #http://finzi.psych.upenn.edu/library/bsts/html/SuggestBurn.html Suggest the size of an MCMC burn in sample as a proportion of the total run.
               seed= 2125)
## =-=-=-=-= Iteration 0 Tue Feb 09 00:11:10 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 3000 Tue Feb 09 00:11:24 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 6000 Tue Feb 09 00:11:39 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 9000 Tue Feb 09 00:11:54 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 12000 Tue Feb 09 00:12:08 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 15000 Tue Feb 09 00:12:23 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 18000 Tue Feb 09 00:12:37 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 21000 Tue Feb 09 00:12:52 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 24000 Tue Feb 09 00:13:06 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 27000 Tue Feb 09 00:13:20 2021
##  =-=-=-=-=
#,
#               dynamic.regression=T)
plot(model1b, main = "Model 1")

plot(model1b, "components")

impact2b <- CausalImpact(bsts.model = model1b,
                       post.period.response = post_period_response)
plot(impact2b, "original") 

burn1b <- SuggestBurn(0.1, model1b)

plot(impact2b)

#summary(impact2b) 
##summary(impact2b,"report")
#msts(data15a64_rn$hosp_trauma, seasonal.periods=c(7,365.25))
#The response variable (i.e., the first column in data) may contain missing values (NA), but covariates (all other columns in data) may not. If one of your covariates contains missing values, consider imputing (i.e., estimating) the missing values; if this is not feasible, leave the regressor out.

# Model 2
ss2b <- list()
# Local trend, weekly-seasonal #https://qastack.mx/stats/209426/predictions-from-bsts-model-in-r-are-failing-completely - PUSE UN GENERALIZED LOCAL TREND
ss2b <- AddLocalLinearTrend(ss2b, data15a64_rn_causal$cons_trauma) #
# Add weekly seasonal
ss2b <- AddSeasonal(ss2b, data15a64_rn_causal$cons_trauma, nseasons=5, season.duration = 52) #weeks OJO, ESTOS NO SON WEEKS VERDADEROS. PORQUE TENGO MAS DE EUN AÑO
ss2b <- AddSeasonal(ss2b, data15a64_rn_causal$cons_trauma, nseasons = 12, season.duration =4) #years
#ss2 <- AddAutoAr(ss2, y = data15a64_rn_causal$log_hosp_trauma, lags = 1) #NO PUEDO AREGAR AR1 CON POISSON
# For example, to add a day-of-week component to data with daily granularity, use model.args = list(nseasons = 7, season.duration = 1). To add a day-of-week component to data with hourly granularity, set model.args = list(nseasons = 7, season.duration = 24).
model2b <- bsts(data15a64_rn_causal$cons_trauma ~ x1,
               state.specification = ss2b, #A list with elements created by AddLocalLinearTrend, AddSeasonal, and similar functions for adding components of state. See the help page for state.specification.
              # family ="poisson", #A Bayesian Analysis of Time-Series Event Count Data
               niter = clus_iter, 
              # burn = 500, #http://finzi.psych.upenn.edu/library/bsts/html/SuggestBurn.html Suggest the size of an MCMC burn in sample as a proportion of the total run.
               seed= 2125)
## =-=-=-=-= Iteration 0 Tue Feb 09 00:14:21 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 3000 Tue Feb 09 00:14:35 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 6000 Tue Feb 09 00:14:50 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 9000 Tue Feb 09 00:15:04 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 12000 Tue Feb 09 00:15:19 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 15000 Tue Feb 09 00:15:33 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 18000 Tue Feb 09 00:15:48 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 21000 Tue Feb 09 00:16:03 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 24000 Tue Feb 09 00:16:17 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 27000 Tue Feb 09 00:16:32 2021
##  =-=-=-=-=
#,
#               dynamic.regression=T)
plot(model2b, main = "Model 2")

plot(model2b, "components")

impact3b <- CausalImpact(bsts.model = model2b,
                       post.period.response = post_period_response)
plot(impact3b, "original") 

burn2b <- SuggestBurn(0.1, model2b)

plot(impact3b)

#summary(impact3b) 
##summary(impact3b,"report")
#dynamic.regression Whether to include time-varying regression coefficients. In combination with a time-varying local trend or even a time-varying local level, this often leads to overspecification, in which case a static regression is safer. Defaults to FALSE.
#
#Prior standard deviation of the Gaussian random walk of the local level. Expressed in terms of data standard deviations. Defaults to 0.01, a typical choice for well-behaved and stable datasets with low residual volatility after regressing out known predictors (e.g., web searches or sales in high quantities). When in doubt, a safer option is to use 0.1, as validated on synthetic data, although this may sometimes give rise to unrealistically wide prediction intervals.


1.1.4 Local Linear Trend & AR


We also introduced an a sparse AR(1) process to the state distribution.


# Model 1
ssb <- list()
# Local trend, weekly-seasonal #https://qastack.mx/stats/209426/predictions-from-bsts-model-in-r-are-failing-completely - PUSE UN GENERALIZED LOCAL TREND
ssb <- AddLocalLinearTrend(ssb, data15a64_rn_causal$cons_trauma) #AddSemilocalLinearTrend #AddLocalLevel
ssb <- AddAutoAr(ssb,data15a64_rn_causal$cons_trauma)

# Add weekly seasonal
ssb <- AddSeasonal(ssb, data15a64_rn_causal$cons_trauma, nseasons=5, season.duration = 52) #weeks OJO, ESTOS NO SON WEEKS VERDADEROS. PORQUE TENGO MAS DE EUN AÑO
ssb <- AddSeasonal(ssb, data15a64_rn_causal$cons_trauma, nseasons = 12, season.duration =4) #years

# For example, to add a day-of-week component to data with daily granularity, use model.args = list(nseasons = 7, season.duration = 1). To add a day-of-week component to data with hourly granularity, set model.args = list(nseasons = 7, season.duration = 24).
model1bar <- bsts(data15a64_rn_causal$cons_trauma, 
               state.specification = ssb, #A list with elements created by AddLocalLinearTrend, AddSeasonal, and similar functions for adding components of state. See the help page for state.specification.
               #family ="poisson", #A Bayesian Analysis of Time-Series Event Count Data. NO SE PUEDE OCUPAR
               niter = clus_iter, 
               #burn = 200, #http://finzi.psych.upenn.edu/library/bsts/html/SuggestBurn.html Suggest the size of an MCMC burn in sample as a proportion of the total run.
               seed= 2125)
## =-=-=-=-= Iteration 0 Tue Feb 09 00:17:38 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 3000 Tue Feb 09 00:17:56 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 6000 Tue Feb 09 00:18:15 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 9000 Tue Feb 09 00:18:33 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 12000 Tue Feb 09 00:18:52 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 15000 Tue Feb 09 00:19:10 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 18000 Tue Feb 09 00:19:29 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 21000 Tue Feb 09 00:19:47 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 24000 Tue Feb 09 00:20:06 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 27000 Tue Feb 09 00:20:24 2021
##  =-=-=-=-=
#,
#               dynamic.regression=T)
plot(model1bar, main = "Model 1")

plot(model1bar, "components")

impact2bar <- CausalImpact(bsts.model = model1bar,
                       post.period.response = post_period_response)
plot(impact2bar, "original") 

burn1bar <- SuggestBurn(0.1, model1bar)

plot(impact2bar)

#summary(impact2b) 
##summary(impact2b,"report")
#msts(data15a64_rn$hosp_trauma, seasonal.periods=c(7,365.25))
#The response variable (i.e., the first column in data) may contain missing values (NA), but covariates (all other columns in data) may not. If one of your covariates contains missing values, consider imputing (i.e., estimating) the missing values; if this is not feasible, leave the regressor out.

# Model 2
ss2b <- list()
# Local trend, weekly-seasonal #https://qastack.mx/stats/209426/predictions-from-bsts-model-in-r-are-failing-completely - PUSE UN GENERALIZED LOCAL TREND
ss2b <- AddLocalLinearTrend(ss2b, data15a64_rn_causal$cons_trauma) #
ss2b <- AddAutoAr(ss2b,data15a64_rn_causal$cons_trauma)

# Add weekly seasonal
ss2b <- AddSeasonal(ss2b, data15a64_rn_causal$cons_trauma, nseasons=5, season.duration = 52) #weeks OJO, ESTOS NO SON WEEKS VERDADEROS. PORQUE TENGO MAS DE EUN AÑO
ss2b <- AddSeasonal(ss2b, data15a64_rn_causal$cons_trauma, nseasons = 12, season.duration =4) #years
#ss2 <- AddAutoAr(ss2, y = data15a64_rn_causal$log_hosp_trauma, lags = 1) #NO PUEDO AREGAR AR1 CON POISSON
# For example, to add a day-of-week component to data with daily granularity, use model.args = list(nseasons = 7, season.duration = 1). To add a day-of-week component to data with hourly granularity, set model.args = list(nseasons = 7, season.duration = 24).
model2bar <- bsts(data15a64_rn_causal$cons_trauma ~ x1,
               state.specification = ss2b, #A list with elements created by AddLocalLinearTrend, AddSeasonal, and similar functions for adding components of state. See the help page for state.specification.
              # family ="poisson", #A Bayesian Analysis of Time-Series Event Count Data
               niter = clus_iter, 
              # burn = 500, #http://finzi.psych.upenn.edu/library/bsts/html/SuggestBurn.html Suggest the size of an MCMC burn in sample as a proportion of the total run.
               seed= 2125)
## =-=-=-=-= Iteration 0 Tue Feb 09 00:21:36 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 3000 Tue Feb 09 00:21:54 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 6000 Tue Feb 09 00:22:13 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 9000 Tue Feb 09 00:22:31 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 12000 Tue Feb 09 00:22:50 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 15000 Tue Feb 09 00:23:09 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 18000 Tue Feb 09 00:23:27 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 21000 Tue Feb 09 00:23:46 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 24000 Tue Feb 09 00:24:04 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 27000 Tue Feb 09 00:24:23 2021
##  =-=-=-=-=
#,
#               dynamic.regression=T)
plot(model2bar, main = "Model 2")

plot(model2bar, "components")

impact3bar <- CausalImpact(bsts.model = model2b,
                       post.period.response = post_period_response)
plot(impact3bar, "original") 

burn2bar <- SuggestBurn(0.1, model2bar)

plot(impact3bar)

#summary(impact3b) 
##summary(impact3b,"report")
#dynamic.regression Whether to include time-varying regression coefficients. In combination with a time-varying local trend or even a time-varying local level, this often leads to overspecification, in which case a static regression is safer. Defaults to FALSE.
#
#Prior standard deviation of the Gaussian random walk of the local level. Expressed in terms of data standard deviations. Defaults to 0.01, a typical choice for well-behaved and stable datasets with low residual volatility after regressing out known predictors (e.g., web searches or sales in high quantities). When in doubt, a safer option is to use 0.1, as validated on synthetic data, although this may sometimes give rise to unrealistically wide prediction intervals.


1.1.5 Prior Local Level SD=.1 (more restrictive)

We used the term prior.level.sd to .1 instead of the default .01, which had been a typical choice for well-behaved and stable datasets with low residual volatility. We decided to distinguish between models that assumed a Random Walk structure, to local linear trends.


# Model 1
ssc <- list()
# Local trend, weekly-seasonal #https://qastack.mx/stats/209426/predictions-from-bsts-model-in-r-are-failing-completely - PUSE UN GENERALIZED LOCAL TREND
ssc <- AddLocalLevel(ssc, data15a64_rn_causal$cons_trauma) #AddSemilocalLinearTrend #AddLocalLevel
# Add weekly seasonal
ssc <- AddSeasonal(ssc, data15a64_rn_causal$cons_trauma, nseasons=5, season.duration = 52) #weeks OJO, ESTOS NO SON WEEKS VERDADEROS. PORQUE TENGO MAS DE EUN AÑO
ssc <- AddSeasonal(ssc, data15a64_rn_causal$cons_trauma, nseasons = 12, season.duration =4) #years
# For example, to add a day-of-week component to data with daily granularity, use model.args = list(nseasons = 7, season.duration = 1). To add a day-of-week component to data with hourly granularity, set model.args = list(nseasons = 7, season.duration = 24).
model1c <- bsts(data15a64_rn_causal$cons_trauma, 
               state.specification = ssc, #A list with elements created by AddLocalLinearTrend, AddSeasonal, and similar functions for adding components of state. See the help page for state.specification.
               #family ="poisson", #A Bayesian Analysis of Time-Series Event Count Data. NO SE PUEDE OCUPAR
               niter = clus_iter, 
               #burn = 200, #http://finzi.psych.upenn.edu/library/bsts/html/SuggestBurn.html Suggest the size of an MCMC burn in sample as a proportion of the total run.
               seed= 2125)
## =-=-=-=-= Iteration 0 Tue Feb 09 00:25:36 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 3000 Tue Feb 09 00:25:48 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 6000 Tue Feb 09 00:26:01 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 9000 Tue Feb 09 00:26:14 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 12000 Tue Feb 09 00:26:26 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 15000 Tue Feb 09 00:26:39 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 18000 Tue Feb 09 00:26:51 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 21000 Tue Feb 09 00:27:04 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 24000 Tue Feb 09 00:27:16 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 27000 Tue Feb 09 00:27:28 2021
##  =-=-=-=-=
#,
#               dynamic.regression=T)
plot(model1c, main = "Model 1")

plot(model1c, "components")

impact2c <- CausalImpact(bsts.model = model1c,model.args = list(prior.level.sd=.1),
                       post.period.response = post_period_response)
plot(impact2c, "original") 

burn1c <- SuggestBurn(0.1, model1c)

plot(impact2c)

#summary(impact2c) 
##summary(impact2c,"report")
#msts(data15a64_rn$hosp_trauma, seasonal.periods=c(7,365.25))
#The response variable (i.e., the first column in data) may contain missing values (NA), but covariates (all other columns in data) may not. If one of your covariates contains missing values, consider imputing (i.e., estimating) the missing values; if this is not feasible, leave the regressor out.

# Model 2
ss2c <- list()
# Local trend, weekly-seasonal #https://qastack.mx/stats/209426/predictions-from-bsts-model-in-r-are-failing-completely - PUSE UN GENERALIZED LOCAL TREND
ss2c <- AddLocalLevel(ss2c, data15a64_rn_causal$cons_trauma) #
# Add weekly seasonal
ss2c <- AddSeasonal(ss2c, data15a64_rn_causal$cons_trauma, nseasons=5, season.duration = 52) #weeks OJO, ESTOS NO SON WEEKS VERDADEROS. PORQUE TENGO MAS DE EUN AÑO
ss2c <- AddSeasonal(ss2c, data15a64_rn_causal$cons_trauma, nseasons = 12, season.duration =4) #years
#ss2 <- AddAutoAr(ss2, y = data15a64_rn_causal$log_hosp_trauma, lags = 1) #NO PUEDO AREGAR AR1 CON POISSON
# For example, to add a day-of-week component to data with daily granularity, use model.args = list(nseasons = 7, season.duration = 1). To add a day-of-week component to data with hourly granularity, set model.args = list(nseasons = 7, season.duration = 24).
model2c <- bsts(data15a64_rn_causal$cons_trauma ~ x1,
               state.specification = ss2c, #A list with elements created by AddLocalLinearTrend, AddSeasonal, and similar functions for adding components of state. See the help page for state.specification.
              # family ="poisson", #A Bayesian Analysis of Time-Series Event Count Data
               niter = clus_iter, 
              # burn = 500, #http://finzi.psych.upenn.edu/library/bsts/html/SuggestBurn.html Suggest the size of an MCMC burn in sample as a proportion of the total run.
               seed= 2125)
## =-=-=-=-= Iteration 0 Tue Feb 09 00:28:21 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 3000 Tue Feb 09 00:28:34 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 6000 Tue Feb 09 00:28:46 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 9000 Tue Feb 09 00:28:58 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 12000 Tue Feb 09 00:29:10 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 15000 Tue Feb 09 00:29:23 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 18000 Tue Feb 09 00:29:35 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 21000 Tue Feb 09 00:29:47 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 24000 Tue Feb 09 00:29:59 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 27000 Tue Feb 09 00:30:12 2021
##  =-=-=-=-=
#,
#               dynamic.regression=T)
plot(model2c, main = "Model 2")

plot(model2c, "components")

impact3c <- CausalImpact(bsts.model = model2c,model.args = list(prior.level.sd=.1),
                       post.period.response = post_period_response)
plot(impact3c, "original") 

burn2c <- SuggestBurn(0.1, model2c)

plot(impact3c)

#summary(impact3c) 
##summary(impact3c,"report")
#dynamic.regression Whether to include time-varying regression coefficients. In combination with a time-varying local trend or even a time-varying local level, this often leads to overspecification, in which case a static regression is safer. Defaults to FALSE.
#
#Prior standard deviation of the Gaussian random walk of the local level. Expressed in terms of data standard deviations. Defaults to 0.01, a typical choice for well-behaved and stable datasets with low residual volatility after regressing out known predictors (e.g., web searches or sales in high quantities). When in doubt, a safer option is to use 0.1, as validated on synthetic data, although this may sometimes give rise to unrealistically wide prediction intervals.

1.1.6 Studentized Distributed Noise (more restrictive)


The default model assumes a Gaussian noise and a Gaussian random walk. In order to handle outliers, we assumed studentized distributed noise.


# Model 1
ssd <- list()
# Local trend, weekly-seasonal #https://qastack.mx/stats/209426/predictions-from-bsts-model-in-r-are-failing-completely - PUSE UN GENERALIZED LOCAL TREND
ssd <- AddLocalLevel(ssd, data15a64_rn_causal$cons_trauma) #AddSemilocalLinearTrend #AddLocalLevel
# Add weekly seasonal
ssd <- AddSeasonal(ssd, data15a64_rn_causal$cons_trauma,nseasons=5, season.duration = 52) #weeks OJO, ESTOS NO SON WEEKS VERDADEROS. PORQUE TENGO MAS DE EUN AÑO
ssd <- AddSeasonal(ssd, data15a64_rn_causal$cons_trauma, nseasons = 12, season.duration =4) #years
# For example, to add a day-of-week component to data with daily granularity, use model.args = list(nseasons = 7, season.duration = 1). To add a day-of-week component to data with hourly granularity, set model.args = list(nseasons = 7, season.duration = 24).
model1d1 <- bsts(data15a64_rn_causal$cons_trauma, 
               state.specification = ssd, #A list with elements created by AddLocalLinearTrend, AddSeasonal, and similar functions for adding components of state. See the help page for state.specification.
               family ="student", #A Bayesian Analysis of Time-Series Event Count Data. POISSON NO SE PUEDE OCUPAR
               niter = clus_iter, 
               #burn = 200, #http://finzi.psych.upenn.edu/library/bsts/html/SuggestBurn.html Suggest the size of an MCMC burn in sample as a proportion of the total run.
               seed= 2125)
## =-=-=-=-= Iteration 0 Tue Feb 09 00:45:43 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 3000 Tue Feb 09 00:45:58 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 6000 Tue Feb 09 00:46:13 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 9000 Tue Feb 09 00:46:28 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 12000 Tue Feb 09 00:46:43 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 15000 Tue Feb 09 00:46:58 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 18000 Tue Feb 09 00:47:14 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 21000 Tue Feb 09 00:47:29 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 24000 Tue Feb 09 00:47:44 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 27000 Tue Feb 09 00:47:59 2021
##  =-=-=-=-=
#,
#               dynamic.regression=T)
plot(model1d1, main = "Model 1")

plot(model1d1, "components")

impact2d1 <- CausalImpact(bsts.model = model1d1,
                       post.period.response = post_period_response)
plot(impact2d1, "original") 

burn1d1 <- SuggestBurn(0.1, model1d1)

plot(impact2d1)

#summary(impact2d1) 
##summary(impact2d1,"report")
#msts(data15a64_rn$hosp_trauma, seasonal.periods=c(7,365.25))
#The response variable (i.e., the first column in data) may contain missing values (NA), but covariates (all other columns in data) may not. If one of your covariates contains missing values, consider imputing (i.e., estimating) the missing values; if this is not feasible, leave the regressor out.

# Model 2
ss2d <- list()
# Local trend, weekly-seasonal #https://qastack.mx/stats/209426/predictions-from-bsts-model-in-r-are-failing-completely - PUSE UN GENERALIZED LOCAL TREND
ss2d <- AddLocalLevel(ss2d, data15a64_rn_causal$cons_trauma) #
# Add weekly seasonal
ss2d <- AddSeasonal(ss2d, data15a64_rn_causal$cons_trauma,nseasons=5, season.duration = 52) #weeks OJO, ESTOS NO SON WEEKS VERDADEROS. PORQUE TENGO MAS DE EUN AÑO
ss2d <- AddSeasonal(ss2d, data15a64_rn_causal$cons_trauma, nseasons = 12, season.duration =4) #years
#ss2 <- AddAutoAr(ss2, y = data15a64_rn_causal$log_hosp_trauma, lags = 1) #NO PUEDO AREGAR AR1 CON POISSON
# For example, to add a day-of-week component to data with daily granularity, use model.args = list(nseasons = 7, season.duration = 1). To add a day-of-week component to data with hourly granularity, set model.args = list(nseasons = 7, season.duration = 24).
model2d1 <- bsts(data15a64_rn_causal$cons_trauma ~ x1,
               state.specification = ss2d, #A list with elements created by AddLocalLinearTrend, AddSeasonal, and similar functions for adding components of state. See the help page for state.specification.
               family ="student", #A Bayesian Analysis of Time-Series Event Count Data
               niter = clus_iter, 
              # burn = 500, #http://finzi.psych.upenn.edu/library/bsts/html/SuggestBurn.html Suggest the size of an MCMC burn in sample as a proportion of the total run.
               seed= 2125)
## =-=-=-=-= Iteration 0 Tue Feb 09 00:48:54 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 3000 Tue Feb 09 00:49:10 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 6000 Tue Feb 09 00:49:25 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 9000 Tue Feb 09 00:49:40 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 12000 Tue Feb 09 00:49:56 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 15000 Tue Feb 09 00:50:11 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 18000 Tue Feb 09 00:50:26 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 21000 Tue Feb 09 00:50:41 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 24000 Tue Feb 09 00:50:56 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 27000 Tue Feb 09 00:51:11 2021
##  =-=-=-=-=
#,
#               dynamic.regression=T)
plot(model2d1, main = "Model 2")

plot(model2d1, "components")

impact3d1 <- CausalImpact(bsts.model = model2d1,
                       post.period.response = post_period_response)
plot(impact3d1, "original") 

burn2d1 <- SuggestBurn(0.1, model2d1)

plot(impact3d1)

#summary(impact3d1) 
##summary(impact3d1,"report")
#dynamic.regression Whether to include time-varying regression coefficients. In combination with a time-varying local trend or even a time-varying local level, this often leads to overspecification, in which case a static regression is safer. Defaults to FALSE.
#
#Prior standard deviation of the Gaussian random walk of the local level. Expressed in terms of data standard deviations. Defaults to 0.01, a typical choice for well-behaved and stable datasets with low residual volatility after regressing out known predictors (e.g., web searches or sales in high quantities). When in doubt, a safer option is to use 0.1, as validated on synthetic data, although this may sometimes give rise to unrealistically wide prediction intervals.

1.1.6.1 Studentized Distributed Noise w/ Priors SD=.1 (more restrictive)


For the following models, we used a prior SD of .1


# Model 1
ssd <- list()
# Local trend, weekly-seasonal #https://qastack.mx/stats/209426/predictions-from-bsts-model-in-r-are-failing-completely - PUSE UN GENERALIZED LOCAL TREND
ssd <- AddLocalLevel(ssd, data15a64_rn_causal$cons_trauma) #AddSemilocalLinearTrend #AddLocalLevel AddLocalLinearTrend
# Add weekly seasonal
ssd <- AddSeasonal(ssd, data15a64_rn_causal$cons_trauma, nseasons=5, season.duration = 52) #weeks OJO, ESTOS NO SON WEEKS VERDADEROS. PORQUE TENGO MAS DE EUN AÑO
ssd <- AddSeasonal(ssd, data15a64_rn_causal$cons_trauma, nseasons = 12, season.duration =4) #years
# For example, to add a day-of-week component to data with daily granularity, use model.args = list(nseasons = 7, season.duration = 1). To add a day-of-week component to data with hourly granularity, set model.args = list(nseasons = 7, season.duration = 24).
model1d <- bsts(data15a64_rn_causal$cons_trauma, 
               state.specification = ssd, #A list with elements created by AddLocalLinearTrend, AddSeasonal, and similar functions for adding components of state. See the help page for state.specification.
               family ="student", #A Bayesian Analysis of Time-Series Event Count Data. POISSON NO SE PUEDE OCUPAR
               niter = clus_iter, 
               #burn = 200, #http://finzi.psych.upenn.edu/library/bsts/html/SuggestBurn.html Suggest the size of an MCMC burn in sample as a proportion of the total run.
               seed= 2125)
## =-=-=-=-= Iteration 0 Tue Feb 09 01:08:48 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 3000 Tue Feb 09 01:09:03 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 6000 Tue Feb 09 01:09:17 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 9000 Tue Feb 09 01:09:32 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 12000 Tue Feb 09 01:09:47 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 15000 Tue Feb 09 01:10:02 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 18000 Tue Feb 09 01:10:17 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 21000 Tue Feb 09 01:10:32 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 24000 Tue Feb 09 01:10:46 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 27000 Tue Feb 09 01:11:01 2021
##  =-=-=-=-=
#,
#               dynamic.regression=T)
plot(model1d, main = "Model 1")

plot(model1d, "components")

impact2d <- CausalImpact(bsts.model = model1d,model.args = list(prior.level.sd=.1),
                       post.period.response = post_period_response)
plot(impact2d, "original") 

burn1d <- SuggestBurn(0.1, model1d)

plot(impact2d)

#summary(impact2d) 
##summary(impact2d,"report")
#msts(data15a64_rn$hosp_trauma, seasonal.periods=c(7,365.25))
#The response variable (i.e., the first column in data) may contain missing values (NA), but covariates (all other columns in data) may not. If one of your covariates contains missing values, consider imputing (i.e., estimating) the missing values; if this is not feasible, leave the regressor out.

# Model 2
ss2d <- list()
# Local trend, weekly-seasonal #https://qastack.mx/stats/209426/predictions-from-bsts-model-in-r-are-failing-completely - PUSE UN GENERALIZED LOCAL TREND
ss2d <- AddLocalLevel(ss2d, data15a64_rn_causal$cons_trauma) #
# Add weekly seasonal
ss2d <- AddSeasonal(ss2d, data15a64_rn_causal$cons_trauma, nseasons=5, season.duration = 52) #weeks OJO, ESTOS NO SON WEEKS VERDADEROS. PORQUE TENGO MAS DE EUN AÑO
ss2d <- AddSeasonal(ss2d, data15a64_rn_causal$cons_trauma, nseasons = 12, season.duration =4) #years
#ss2 <- AddAutoAr(ss2, y = data15a64_rn_causal$log_hosp_trauma, lags = 1) #NO PUEDO AREGAR AR1 CON POISSON
# For example, to add a day-of-week component to data with daily granularity, use model.args = list(nseasons = 7, season.duration = 1). To add a day-of-week component to data with hourly granularity, set model.args = list(nseasons = 7, season.duration = 24).
model2d <- bsts(data15a64_rn_causal$cons_trauma ~ x1,
               state.specification = ss2d, #A list with elements created by AddLocalLinearTrend, AddSeasonal, and similar functions for adding components of state. See the help page for state.specification.
               family ="student", #A Bayesian Analysis of Time-Series Event Count Data
               niter = clus_iter, 
              # burn = 500, #http://finzi.psych.upenn.edu/library/bsts/html/SuggestBurn.html Suggest the size of an MCMC burn in sample as a proportion of the total run.
               seed= 2125)
## =-=-=-=-= Iteration 0 Tue Feb 09 01:11:55 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 3000 Tue Feb 09 01:12:10 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 6000 Tue Feb 09 01:12:25 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 9000 Tue Feb 09 01:12:40 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 12000 Tue Feb 09 01:12:56 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 15000 Tue Feb 09 01:13:11 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 18000 Tue Feb 09 01:13:26 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 21000 Tue Feb 09 01:13:42 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 24000 Tue Feb 09 01:13:57 2021
##  =-=-=-=-=
## =-=-=-=-= Iteration 27000 Tue Feb 09 01:14:12 2021
##  =-=-=-=-=
#,
#               dynamic.regression=T)
plot(model2d, main = "Model 2")

plot(model2d, "components")

impact3d <- CausalImpact(bsts.model = model2d,model.args = list(prior.level.sd=.1),
                       post.period.response = post_period_response)
plot(impact3d, "original") 

burn2d <- SuggestBurn(0.1, model2d)

plot(impact3d)

#summary(impact3d) 
##summary(impact3d,"report")
#dynamic.regression Whether to include time-varying regression coefficients. In combination with a time-varying local trend or even a time-varying local level, this often leads to overspecification, in which case a static regression is safer. Defaults to FALSE.
#
#Prior standard deviation of the Gaussian random walk of the local level. Expressed in terms of data standard deviations. Defaults to 0.01, a typical choice for well-behaved and stable datasets with low residual volatility after regressing out known predictors (e.g., web searches or sales in high quantities). When in doubt, a safer option is to use 0.1, as validated on synthetic data, although this may sometimes give rise to unrealistically wide prediction intervals.

1.2 Comparison Between Models


We compared the models selected in terms of cumulative absolute error.


comp_mod<-CompareBstsModels(list("No control variables(a)" = model1,
                                "Control variables(b)" = model2,
                                "a, LocalLinearTrend" = model1b,
                                "b, LocalLinearTrend" = model2b,
                                "a, LocalLinearTrend & AR" = model1bar,
                                "b, LocalLinearTrend & AR" = model2bar,
                                "a, Prior sd=.1" = model1c,
                                "b, Prior sd=.1" = model2c,
                                "a, Prior sd=.1,LocalLinearTrend" = model1c3,
                                "b, Prior sd=.1,LocalLinearTrend" = model2c3,
                                "a, Prior sd=.1,LocalLinearTrend & AR" = model1c3ar,
                                "b, Prior sd=.1,LocalLinearTrend & AR" = model2c3ar,
                                "a, Student Dist, Prior sd=.01" = model1d1, #Student SD=.1 Local
                                "b, Student Dist, Prior sd=.01" = model2d1, #Student SD=.1 Local
                                "a, Student Dist, Prior sd=.01, LocalLinearTrend" = model1d2,
                                "b, Student Dist, Prior sd=.01, LocalLinearTrend" = model2d2,
                                "a, Student Dist, Prior sd=.01, LocalLinearTrend & AR" = model1d2ar,
                                "b, Student Dist, Prior sd=.01, LocalLinearTrend & AR" = model2d2ar,
                                "a, Student Dist, Prior sd=.1" = model1d,
                                "b, Student Dist, Prior sd=.1" = model2d,
                                "a, Student Dist, Prior sd=.1, LocalLinearTrend" = model1d3,
                                "b, Student Dist, Prior sd=.1, LocalLinearTrend" = model2d3,
                                "a, Student Dist, Prior sd=.1, LocalLinearTrend & AR" = model1d3ar,
                                "b, Student Dist, Prior sd=.1, LocalLinearTrend & AR" = model2d3ar
                       ),
                  colors = c("black", "red", "blue","purple","gray","darkslategrey","darkslateblue","darkorange3","forestgreen",
                             "yellow","pink", "brown","cyan"))
scale2 <- function(x, na.rm = FALSE) (x - mean(x, na.rm = na.rm)) / sd(x, na.rm)

comp_mod_df<-
data.frame(cbind(Models=c("No control variables(a)",
                          "Control variables(b)",
                          "a, LocalLinearTrend",
                          "b, LocalLinearTrend",
                          "a, LocalLinearTrend & AR",
                          "b, LocalLinearTrend & AR",                          
                          "a, Prior sd=.1",
                          "b, Prior sd=.1",
                          "a, Prior sd=.1, LocalLinearTrend",
                          "b, Prior sd=.1, LocalLinearTrend",
                          "a, Prior sd=.1, LocalLinearTrend & AR",
                          "b, Prior sd=.1, LocalLinearTrend & AR",
                          "a, Student Dist, Prior sd=.01",
                          "b, Student Dist, Prior sd=.01",
                          "a, Student Dist, Prior sd=.01, LocalLinearTrend",
                          "b, Student Dist, Prior sd=.01, LocalLinearTrend",
                          "a, Student Dist, Prior sd=.01, LocalLinearTrend & AR",
                          "b, Student Dist, Prior sd=.01, LocalLinearTrend & AR",
                          "a, Student Dist, Prior sd=.1",
                          "b, Student Dist, Prior sd=.1",
                          "a, Student Dist, Prior sd=.1, LocalLinearTrend",
                          "b, Student Dist, Prior sd=.1, LocalLinearTrend",
                          "a, Student Dist, Prior sd=.1, LocalLinearTrend & AR",
                          "b, Student Dist, Prior sd=.1, LocalLinearTrend & AR"),
                 comp_mod)) %>%  #16
  melt(id=1)%>%
  dplyr::rename("Time"="variable") %>% 
  dplyr::rename("Cumulative Absolute Error"="value") %>% 
  dplyr::mutate(Time=as.numeric(sub('V', '', Time))) %>% 
  dplyr::mutate(`Cumulative Absolute Error`=as.numeric(`Cumulative Absolute Error`)) %>% 
  dplyr::mutate(text=paste0("CAFE= ",sprintf("%4.0f",`Cumulative Absolute Error`),"\n","Time= ",Time,"\n",
                            Models))

ggplotcomp_models<-    
ggplot(comp_mod_df)+
  #geom_point(color="white")+
    geom_line(size=.75, aes(x = Time, y = `Cumulative Absolute Error`, color=Models,group=Models, text=text))+
    theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  scale_fill_brewer(type="seq", palette="Greys")+
  theme_sjplot2()+
  theme(axis.text.y=element_blank(),
                    axis.text.x =  element_text(angle = 0, hjust = 1),
                    axis.ticks.y=element_blank(),                    
                    panel.grid.major.y=element_blank(),
        #panel.grid.major.x=element_blank(),
        panel.grid.minor.x=element_blank())+
  theme(legend.position='none')+
  scale_x_continuous(breaks=seq(from =1,to =max(comp_mod_df$Time)+4,by=12))+
  labs(x="Time (in weeks)", y= "Cumulative Absolute Forecast Error (CAFE)")

ggplotly(ggplotcomp_models,tooltip="text")

Figure 9. Comparison of BSTS models (cum. error) 

comp_mod_df_sum<-
comp_mod_df %>% group_by(Models) %>% summarise(mean=mean(abs(`Cumulative Absolute Error`),na.rm=T),median=median(abs(`Cumulative Absolute Error`),na.rm=T), p25=quantile(abs(`Cumulative Absolute Error`),.25,na.rm=T),p75=quantile(abs(`Cumulative Absolute Error`),.75,na.rm=T), max=max(`Cumulative Absolute Error`,na.rm=T)) %>%
  dplyr::mutate(mean=scale2(mean),median=scale2(median),p25=scale2(p25),p75=scale2(p75)) %>% 
  rowwise() %>%
  mutate(mean_tot = mean(c(mean,median,p25,p75))) %>% 
  dplyr::arrange(mean_tot)
  #geom_line(aes(x=Models,y=median),color="blue")+
  #geom_line(aes(x=Models,y=p25),color="green")


#impact "No counterfactual approach, No seasonal components" "Counterfactual Approach, NoSeasonal Components (impact2)" (b= SemilocalLinearTrend) (c=Gaussian, prior sd=.01) (d1= Student, Prior sd .01) (d= Student, Prior sd = .1)
#impact2d$summary[,"p"]


As seen in the Figure above, the errors in the different models were very similar. However, some exhibited lower cumulative errors.


summary<-cbind(Models=c("No control variables(a)",
                          "Control variables(b)",
                          "a, LocalLinearTrend",
                          "b, LocalLinearTrend",
                          "a, LocalLinearTrend & AR",
                          "b, LocalLinearTrend & AR",                          
                          "a, Prior sd=.1",
                          "b, Prior sd=.1",
                          "a, Prior sd=.1, LocalLinearTrend",
                          "b, Prior sd=.1, LocalLinearTrend",
                          "a, Prior sd=.1, LocalLinearTrend & AR",
                          "b, Prior sd=.1, LocalLinearTrend & AR",
                          "a, Student Dist, Prior sd=.01",
                          "b, Student Dist, Prior sd=.01",
                          "a, Student Dist, Prior sd=.01, LocalLinearTrend",
                          "b, Student Dist, Prior sd=.01, LocalLinearTrend",
                          "a, Student Dist, Prior sd=.01, LocalLinearTrend & AR",
                          "b, Student Dist, Prior sd=.01, LocalLinearTrend & AR",
                          "a, Student Dist, Prior sd=.1",
                          "b, Student Dist, Prior sd=.1",
                          "a, Student Dist, Prior sd=.1, LocalLinearTrend",
                          "b, Student Dist, Prior sd=.1, LocalLinearTrend",
                          "a, Student Dist, Prior sd=.1, LocalLinearTrend & AR",
                          "b, Student Dist, Prior sd=.1, LocalLinearTrend & AR"),
               name=c("impact2","impact3","impact2b","impact3b","impact2bar","impact3bar","impact2c","impact3c","impact2c3","impact3c3","impact2c3ar","impact3c3ar","impact2d1","impact3d1","impact2d2","impact2d2ar","impact3d2","impact3d2ar","impact2d","impact3d","impact2d3","impact3d3","impact2d3ar","impact3d3ar"),
  `_`= rbind(    
                        data.table(impact2$summary,keep.rownames = T)[1,],
                        data.table(impact3$summary,keep.rownames = T)[1,],
                        data.table(impact2b$summary,keep.rownames = T)[1,],
                        data.table(impact3b$summary,keep.rownames = T)[1,],
                        data.table(impact2bar$summary,keep.rownames = T)[1,],
                        data.table(impact3bar$summary,keep.rownames = T)[1,],
                        data.table(impact2c$summary,keep.rownames = T)[1,],
                        data.table(impact3c$summary,keep.rownames = T)[1,],
                        data.table(impact2c3$summary,keep.rownames = T)[1,],
                        data.table(impact3c3$summary,keep.rownames = T)[1,],
                        data.table(impact2c3ar$summary,keep.rownames = T)[1,],
                        data.table(impact3c3ar$summary,keep.rownames = T)[1,],
                        data.table(impact2d1$summary,keep.rownames = T)[1,],
                        data.table(impact3d1$summary,keep.rownames = T)[1,],
                        data.table(impact2d2$summary,keep.rownames = T)[1,],
                        data.table(impact2d2ar$summary,keep.rownames = T)[1,],
                        data.table(impact3d2$summary,keep.rownames = T)[1,],
                        data.table(impact3d2ar$summary,keep.rownames = T)[1,],
                        data.table(impact2d$summary,keep.rownames = T)[1,],
                        data.table(impact3d$summary,keep.rownames = T)[1,],
                        data.table(impact2d3$summary,keep.rownames = T)[1,],
                        data.table(impact3d3$summary,keep.rownames = T)[1,],
                        data.table(impact2d3ar$summary,keep.rownames = T)[1,],
                        data.table(impact3d3ar$summary,keep.rownames = T)[1,]
             ))%>% data.table(.,keep.rownames = F)

#comp_mod_df_sum[,c("Models","name","mean_tot")]
  summary%>%
    rename_all(~sub('_.', '', .x))%>%
    dplyr::mutate_at(4:(ncol(.)-1),~round(as.numeric(.),2)) %>% 
    dplyr::mutate_at(ncol(.),~round(as.numeric(.),4)) %>% 
   # dplyr::rename_at(vars(starts_with('X_')),~
  knitr::kable(format = "html", format.args = list(decimal.mark = ".", big.mark = ","),
               caption="Table 1. Summary of the effects of the different models",
                 align =c('l',rep('c', 101)))%>%
   kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 8) %>%
  kableExtra::scroll_box(width = "100%", height = "375px")
Table 1. Summary of the effects of the different models
Models name rn Actual Pred Pred.lower Pred.upper Pred.sd AbsEffect AbsEffect.lower AbsEffect.upper AbsEffect.sd RelEffect RelEffect.lower RelEffect.upper RelEffect.sd alpha p
No control variables(a) impact2 Average 806.9 988.86 709.91 1,260.71 140.22 -181.96 -453.81 96.99 140.22 -0.18 -0.46 0.10 0.14 0.05 0.0953
Control variables(b) impact3 Average 806.9 928.32 662.27 1,191.87 135.42 -121.42 -384.97 144.63 135.42 -0.13 -0.41 0.16 0.15 0.05 0.1819
a, LocalLinearTrend impact2b Average 806.9 996.23 703.66 1,285.93 148.39 -189.33 -479.03 103.24 148.39 -0.19 -0.48 0.10 0.15 0.05 0.0980
b, LocalLinearTrend impact3b Average 806.9 931.50 647.73 1,213.54 143.39 -124.60 -406.64 159.17 143.39 -0.13 -0.44 0.17 0.15 0.05 0.1867
a, LocalLinearTrend & AR impact2bar Average 806.9 998.54 776.24 1,212.98 111.27 -191.64 -406.08 30.66 111.27 -0.19 -0.41 0.03 0.11 0.05 0.0452
b, LocalLinearTrend & AR impact3bar Average 806.9 931.50 648.68 1,212.02 143.23 -124.60 -405.12 158.22 143.23 -0.13 -0.43 0.17 0.15 0.05 0.1868
a, Prior sd=.1 impact2c Average 806.9 988.86 710.29 1,261.31 140.01 -181.96 -454.41 96.61 140.01 -0.18 -0.46 0.10 0.14 0.05 0.0942
b, Prior sd=.1 impact3c Average 806.9 928.32 660.62 1,192.18 135.63 -121.42 -385.28 146.28 135.63 -0.13 -0.42 0.16 0.15 0.05 0.1816
a, Prior sd=.1, LocalLinearTrend impact2c3 Average 806.9 996.23 700.64 1,285.21 148.42 -189.33 -478.31 106.26 148.42 -0.19 -0.48 0.11 0.15 0.05 0.0977
b, Prior sd=.1, LocalLinearTrend impact3c3 Average 806.9 931.50 646.40 1,211.54 143.17 -124.60 -404.64 160.50 143.17 -0.13 -0.43 0.17 0.15 0.05 0.1865
a, Prior sd=.1, LocalLinearTrend & AR impact2c3ar Average 806.9 998.54 776.28 1,214.40 111.31 -191.64 -407.50 30.62 111.31 -0.19 -0.41 0.03 0.11 0.05 0.0449
b, Prior sd=.1, LocalLinearTrend & AR impact3c3ar Average 806.9 936.20 721.57 1,143.51 107.37 -129.30 -336.61 85.33 107.37 -0.14 -0.36 0.09 0.11 0.05 0.1142
a, Student Dist, Prior sd=.01 impact2d1 Average 806.9 951.48 691.16 1,211.89 132.07 -144.58 -404.99 115.74 132.07 -0.15 -0.43 0.12 0.14 0.05 0.1341
b, Student Dist, Prior sd=.01 impact3d1 Average 806.9 905.90 645.57 1,164.60 131.29 -99.00 -357.70 161.33 131.29 -0.11 -0.39 0.18 0.14 0.05 0.2206
a, Student Dist, Prior sd=.01, LocalLinearTrend impact2d2 Average 806.9 953.94 683.28 1,227.09 139.35 -147.04 -420.19 123.62 139.35 -0.15 -0.44 0.13 0.15 0.05 0.1468
b, Student Dist, Prior sd=.01, LocalLinearTrend impact2d2ar Average 806.9 953.58 758.65 1,146.67 99.17 -146.68 -339.77 48.25 99.17 -0.15 -0.36 0.05 0.10 0.05 0.0689
a, Student Dist, Prior sd=.01, LocalLinearTrend & AR impact3d2 Average 806.9 909.20 635.63 1,182.45 140.24 -102.30 -375.55 171.27 140.24 -0.11 -0.41 0.19 0.15 0.05 0.2318
b, Student Dist, Prior sd=.01, LocalLinearTrend & AR impact3d2ar Average 806.9 912.02 714.33 1,110.29 100.27 -105.12 -303.39 92.57 100.27 -0.12 -0.33 0.10 0.11 0.05 0.1434
a, Student Dist, Prior sd=.1 impact2d Average 806.9 951.48 691.48 1,212.18 132.04 -144.58 -405.28 115.42 132.04 -0.15 -0.43 0.12 0.14 0.05 0.1326
b, Student Dist, Prior sd=.1 impact3d Average 806.9 905.90 646.30 1,164.99 131.32 -99.00 -358.09 160.60 131.32 -0.11 -0.40 0.18 0.14 0.05 0.2208
a, Student Dist, Prior sd=.1, LocalLinearTrend impact2d3 Average 806.9 953.94 683.02 1,227.46 139.35 -147.04 -420.56 123.88 139.35 -0.15 -0.44 0.13 0.15 0.05 0.1460
b, Student Dist, Prior sd=.1, LocalLinearTrend impact3d3 Average 806.9 909.20 635.76 1,182.81 140.23 -102.30 -375.91 171.14 140.23 -0.11 -0.41 0.19 0.15 0.05 0.2298
a, Student Dist, Prior sd=.1, LocalLinearTrend & AR impact2d3ar Average 806.9 953.58 758.30 1,147.70 99.25 -146.68 -340.80 48.60 99.25 -0.15 -0.36 0.05 0.10 0.05 0.0689
b, Student Dist, Prior sd=.1, LocalLinearTrend & AR impact3d3ar Average 806.9 912.02 714.62 1,109.61 100.21 -105.12 -302.71 92.28 100.21 -0.12 -0.33 0.10 0.11 0.05 0.1428 


comp_mod_df_sum<-
  comp_mod_df_sum %>% 
    dplyr::filter(grepl("^a,", Models)|grepl("\\(a,", Models)) %>% 
    dplyr::left_join(summary[,c("Models","name")],by="Models")

comp_mod_df_sum %>% 
  dplyr::filter(grepl("^a,", Models)) %>% 
  dplyr::select( Models, mean, median, p25, p75, mean_tot, max) %>% 
    dplyr::mutate_at(2:(ncol(.)-1),~round(as.numeric(.),2)) %>% 
      dplyr::mutate_at(ncol(.),~round(as.numeric(.),0)) %>% 
    #dplyr::arrange(Models,rn) %>% 
   # dplyr::rename_at(vars(starts_with('X_')),~rstudio
  knitr::kable(format = "html", format.args = list(decimal.mark = ".", big.mark = ","),
               caption="Table 2. Standardized Comparison of Estimated Models",
               col.names=c("Models", "Mean (Std. z scale)","Median (z scale)", "Q1 (z scale)", "Q3 (z scale)",
                           "Total Average (z scale)","Maximum CAFE"),
                 align =c('l',rep('c', 101)))%>%
   kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 8) %>%
  kableExtra::scroll_box(width = "100%", height = "375px")
Table 2. Standardized Comparison of Estimated Models
Models Mean (Std. z scale) Median (z scale) Q1 (z scale) Q3 (z scale) Total Average (z scale) Maximum CAFE
a, Student Dist, Prior sd=.01 -0.84 -0.81 -0.45 -0.97 -0.77 16,348
a, Student Dist, Prior sd=.1 -0.84 -0.81 -0.45 -0.97 -0.77 16,348
a, Student Dist, Prior sd=.01, LocalLinearTrend & AR -0.59 -0.62 0.19 -0.88 -0.47 16,263
a, Student Dist, Prior sd=.1, LocalLinearTrend & AR -0.59 -0.62 0.19 -0.88 -0.47 16,263
a, Student Dist, Prior sd=.01, LocalLinearTrend -0.36 -0.44 0.15 -0.58 -0.31 16,694
a, Student Dist, Prior sd=.1, LocalLinearTrend -0.36 -0.44 0.15 -0.58 -0.31 16,694
a, Prior sd=.1 0.16 0.37 -0.25 0.22 0.12 17,643
a, LocalLinearTrend 1.56 1.60 1.57 1.31 1.51 18,540
a, Prior sd=.1, LocalLinearTrend 1.56 1.60 1.57 1.31 1.51 18,540
a, LocalLinearTrend & AR 1.71 1.74 2.01 1.46 1.73 18,322
a, Prior sd=.1, LocalLinearTrend & AR 1.71 1.74 2.01 1.46 1.73 18,322

1.3 Selected Model


For the meantime, we selected the model with the lowest cumulative errors.


#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
invisible(c("'0b, Student Dist, Prior sd=.01' = model2d1, #Student SD=.1 Local           'b, Prior sd=.1'"))

cat("#### ",as.character(comp_mod_df_sum[1,c("Models")]),"\n")
## ####  a, Student Dist, Prior sd=.01
summary(get(as.character(comp_mod_df_sum[1,c("name")]))) 
## Posterior inference {CausalImpact}
## 
##                          Average       Cumulative   
## Actual                   807           8069         
## Prediction (s.d.)        951 (132)     9515 (1321)  
## 95% CI                   [691, 1212]   [6912, 12119]
##                                                     
## Absolute effect (s.d.)   -145 (132)    -1446 (1321) 
## 95% CI                   [-405, 116]   [-4050, 1157]
##                                                     
## Relative effect (s.d.)   -15% (14%)    -15% (14%)   
## 95% CI                   [-43%, 12%]   [-43%, 12%]  
## 
## Posterior tail-area probability p:   0.13414
## Posterior prob. of a causal effect:  87%
## 
## For more details, type: summary(impact, "report")
cat("#### ",as.character(comp_mod_df_sum[2,c("Models")]),"\n")
## ####  a, Student Dist, Prior sd=.1
summary(get(as.character(comp_mod_df_sum[2,c("name")])))
## Posterior inference {CausalImpact}
## 
##                          Average       Cumulative   
## Actual                   807           8069         
## Prediction (s.d.)        951 (132)     9515 (1320)  
## 95% CI                   [691, 1212]   [6915, 12122]
##                                                     
## Absolute effect (s.d.)   -145 (132)    -1446 (1320) 
## 95% CI                   [-405, 115]   [-4053, 1154]
##                                                     
## Relative effect (s.d.)   -15% (14%)    -15% (14%)   
## 95% CI                   [-43%, 12%]   [-43%, 12%]  
## 
## Posterior tail-area probability p:   0.13262
## Posterior prob. of a causal effect:  87%
## 
## For more details, type: summary(impact, "report")
#The “Average” column represents the average across time during the post-intervention period, while the “Cumulative” column presents the total sum of the time points. In particular, looking at the Average Absolute effect, we see that it is estimated to be 22, with a 95% posterior interval between -38 to -5. Since the interval excludes 0, we can conclude that the intervention of the Public Health Emergency in Los Angeles had a causal impact on the PM2.5 air quality with certain assumptions.



2 References

  • Cox, L. (2015) Quantifying and Reducing Uncertainty about Causality in Improving Public Health and Safety. In: Ghanem R., Higdon D., Owhadi H. (eds) Handbook of Uncertainty Quantification. Springer, Cham. https://doi.org/10.1007/978-3-319-11259-6_71-1

  • Brodersen, KH., Gallusser, F., Koehler, J., Remy, N., Scott, SL. (2015) Inferring causal impact using Bayesian structural time-series models. Ann Appl Stat 9:247-274